/***********************************************************************
/
/  INITIALIZE A COSMOLOGY SIMULATION
/
/  written by: Greg Bryan
/  date:       May, 1995
/  modified1:  Robert Harkness
/  date:       July 2003
/
/  PURPOSE:  Initialize for cosmology simulations.  Reads in a number
/      of initial grids.  If more than one, then they are all numbered.
/      We currently assume that all are subgrids except the first.
/
/  RETURNS: SUCCESS or FAIL
/
************************************************************************/
 
// This routine intializes a new simulation based on the parameter file
 
 
#ifdef USE_MPI
#include "mpi.h"
#endif /* USE_MPI */
 
#include <string.h>
#include <stdio.h>
#include <math.h>
 
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
#include "Fluxes.h"
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
#include "CosmologyParameters.h"
#include "fortran.def"
#include "CommunicationUtilities.h"

// Function prototypes
 
void WriteListOfFloats(FILE *fptr, int N, float floats[]);
void WriteListOfFloats(FILE *fptr, int N, FLOAT floats[]);
void WriteListOfInts(FILE *fptr, int N, int nums[]);
int CommunicationBroadcastValue(PINT *Value, int BroadcastProcessor);
void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
 
// Cosmology Parameters (that need to be shared)
 
static float CosmologySimulationOmegaBaryonNow       = 1.0;  // standard
static float CosmologySimulationOmegaCDMNow          = 0.0;  // no dark matter
static float CosmologySimulationInitialTemperature   = FLOAT_UNDEFINED;
 
static char *CosmologySimulationDensityName          = NULL;
static char *CosmologySimulationTotalEnergyName      = NULL;
static char *CosmologySimulationGasEnergyName        = NULL;
static char *CosmologySimulationParticlePositionName = NULL;
static char *CosmologySimulationParticleVelocityName = NULL;
static char *CosmologySimulationParticleMassName     = NULL;
static char *CosmologySimulationParticleTypeName     = NULL;
static char *CosmologySimulationVelocityNames[MAX_DIMENSION];
static char *CosmologySimulationParticleVelocityNames[MAX_DIMENSION];
 
static int   CosmologySimulationSubgridsAreStatic    = TRUE;
 
static float CosmologySimulationInitialFractionHII   = 1.2e-5;
static float CosmologySimulationInitialFractionHeII  = 1.0e-14;
static float CosmologySimulationInitialFractionHeIII = 1.0e-17;
static float CosmologySimulationInitialFractionHM    = 2.0e-9;
static float CosmologySimulationInitialFractionH2I   = 2.0e-20;
static float CosmologySimulationInitialFractionH2II  = 3.0e-14;
static int   CosmologySimulationUseMetallicityField  = FALSE;
 
static int CosmologySimulationManuallySetParticleMassRatio = FALSE;
static float CosmologySimulationManualParticleMassRatio = 1.0;

static int   CosmologySimulationCalculatePositions   = FALSE; 

#define MAX_INITIAL_GRIDS 10
 
 
int NestedCosmologySimulationInitialize(FILE *fptr, FILE *Outfptr,
			       		HierarchyEntry &TopGrid, TopGridData &MetaData)
{
  char *DensName = "Density";
  char *TEName   = "TotalEnergy";
  char *GEName   = "GasEnergy";
  char *Vel1Name = "x-velocity";
  char *Vel2Name = "y-velocity";
  char *Vel3Name = "z-velocity";
  char *ElectronName = "Electron_Density";
  char *HIName    = "HI_Density";
  char *HIIName   = "HII_Density";
  char *HeIName   = "HeI_Density";
  char *HeIIName  = "HeII_Density";
  char *HeIIIName = "HeIII_Density";
  char *HMName    = "HM_Density";
  char *H2IName   = "H2I_Density";
  char *H2IIName  = "H2II_Density";
  char *DIName    = "DI_Density";
  char *DIIName   = "DII_Density";
  char *HDIName   = "HDI_Density";
  char *GPotName  = "Grav_Potential";
  char *MetalName = "Metal_Density";
  char *ForbidName = "ForbiddenRefinement";
 
  char *ExtraNames[2] = {"Z_Field1", "Z_Field2"};
 
  // Declarations
 
  char line[MAX_LINE_LENGTH];
  int i, j, dim, gridnum, ret, SubgridsAreStatic, region;
  HierarchyEntry *Subgrid;
 
  char *DensityName = NULL, *TotalEnergyName = NULL, *GasEnergyName = NULL,
    *ParticlePositionName = NULL, *ParticleVelocityName = NULL,
    *ParticleMassName = NULL, *VelocityNames[MAX_DIMENSION],
    *ParticleTypeName = NULL, *ParticleVelocityNames[MAX_DIMENSION];
 
  for (dim = 0; dim < MAX_DIMENSION; dim++) {
    VelocityNames[dim] = NULL;
    ParticleVelocityNames[dim] = NULL;
    CosmologySimulationParticleVelocityNames[dim] = NULL;
  }
 
  // Set default parameters: parameters, names and subgrid info
 
  for (dim = 0; dim < MAX_DIMENSION; dim++)
    CosmologySimulationVelocityNames[dim]       = NULL;
 
  int   CosmologySimulationGridDimension[MAX_INITIAL_GRIDS][MAX_DIMENSION];
  int   CosmologySimulationGridLevel[MAX_INITIAL_GRIDS];
  FLOAT CosmologySimulationGridLeftEdge[MAX_INITIAL_GRIDS][MAX_DIMENSION];
  FLOAT CosmologySimulationGridRightEdge[MAX_INITIAL_GRIDS][MAX_DIMENSION];
 
  for (i = 0; i < MAX_INITIAL_GRIDS; i++)
    CosmologySimulationGridLevel[i] = 1;
 
  for (dim = 0; dim < MetaData.TopGridRank; dim++) {
    CosmologySimulationGridLeftEdge[0][dim] = DomainLeftEdge[dim];
    CosmologySimulationGridRightEdge[0][dim] = DomainRightEdge[dim];
    CosmologySimulationGridDimension[0][dim] = MetaData.TopGridDims[dim];
  }
 
  CosmologySimulationGridLevel[0] = 0;
 
  // Error check
 
  if (!ComovingCoordinates) {
    ENZO_FAIL("ComovingCoordinates must be TRUE!\n");
  }
 
  if (DualEnergyFormalism == FALSE && HydroMethod != Zeus_Hydro)
    fprintf(stderr, "CosmologySimulation: DualEnergyFormalism is off!\n");
  if (!SelfGravity)
    fprintf(stderr, "CosmologySimulation: gravity is off!?!\n");
 
  // Read keyword input from file
 
  char *dummy = new char[MAX_LINE_LENGTH];
 
  dummy[0] = 0;
 
  while (fgets(line, MAX_LINE_LENGTH, fptr) != NULL) {
 
    ret = 0;
 
    ret += sscanf(line, "CosmologySimulationOmegaBaryonNow = %"FSYM,
		  &CosmologySimulationOmegaBaryonNow);
    ret += sscanf(line, "CosmologySimulationOmegaCDMNow = %"FSYM,
		  &CosmologySimulationOmegaCDMNow);
    ret += sscanf(line, "CosmologySimulationInitialTemperature = %"FSYM,
		  &CosmologySimulationInitialTemperature);
 
    if (sscanf(line, "CosmologySimulationDensityName = %s", dummy) == 1)
      CosmologySimulationDensityName = dummy;
    if (sscanf(line, "CosmologySimulationTotalEnergyName = %s", dummy) == 1)
      CosmologySimulationTotalEnergyName = dummy;
    if (sscanf(line, "CosmologySimulationGasEnergyName = %s", dummy) == 1)
      CosmologySimulationGasEnergyName = dummy;
    if (sscanf(line, "CosmologySimulationVelocity1Name = %s", dummy) == 1)
      CosmologySimulationVelocityNames[0] = dummy;
    if (sscanf(line, "CosmologySimulationVelocity2Name = %s", dummy) == 1)
      CosmologySimulationVelocityNames[1] = dummy;
    if (sscanf(line, "CosmologySimulationVelocity3Name = %s", dummy) == 1)
      CosmologySimulationVelocityNames[2] = dummy;
    if (sscanf(line, "CosmologySimulationParticlePositionName = %s", dummy) == 1)
      CosmologySimulationParticlePositionName = dummy;
    if (sscanf(line, "CosmologySimulationParticleVelocityName = %s", dummy) == 1)
      CosmologySimulationParticleVelocityName = dummy;
    if (sscanf(line, "CosmologySimulationParticleMassName = %s", dummy) == 1)
      CosmologySimulationParticleMassName = dummy;
    if (sscanf(line, "CosmologySimulationParticleTypeName = %s", dummy) == 1)
      CosmologySimulationParticleTypeName = dummy;
    if (sscanf(line, "CosmologySimulationParticleVelocity1Name = %s", dummy) == 1)
      CosmologySimulationParticleVelocityNames[0] = dummy;
    if (sscanf(line, "CosmologySimulationParticleVelocity2Name = %s", dummy) == 1)
      CosmologySimulationParticleVelocityNames[1] = dummy;
    if (sscanf(line, "CosmologySimulationParticleVelocity3Name = %s", dummy) == 1)
      CosmologySimulationParticleVelocityNames[2] = dummy;    

    ret += sscanf(line, "CosmologySimulationNumberOfInitialGrids = %"ISYM,
		  &CosmologySimulationNumberOfInitialGrids);
    ret += sscanf(line, "CosmologySimulationSubgridsAreStatic = %"ISYM,
		  &CosmologySimulationSubgridsAreStatic);
 
    if (sscanf(line, "CosmologySimulationGridLeftEdge[%"ISYM"]", &gridnum) > 0)
      ret += sscanf(line, "CosmologySimulationGridLeftEdge[%"ISYM"] = %"PSYM" %"PSYM" %"PSYM,
		    &gridnum, &CosmologySimulationGridLeftEdge[gridnum][0],
		    &CosmologySimulationGridLeftEdge[gridnum][1],
		    &CosmologySimulationGridLeftEdge[gridnum][2]);
    if (sscanf(line, "CosmologySimulationGridRightEdge[%"ISYM"]", &gridnum) > 0)
      ret += sscanf(line, "CosmologySimulationGridRightEdge[%"ISYM"] = %"PSYM" %"PSYM" %"PSYM,
		    &gridnum, &CosmologySimulationGridRightEdge[gridnum][0],
		    &CosmologySimulationGridRightEdge[gridnum][1],
		    &CosmologySimulationGridRightEdge[gridnum][2]);
    if (sscanf(line, "CosmologySimulationGridDimension[%"ISYM"]", &gridnum) > 0)
      ret += sscanf(line, "CosmologySimulationGridDimension[%"ISYM"] = %"ISYM" %"ISYM" %"ISYM,
		    &gridnum, &CosmologySimulationGridDimension[gridnum][0],
		    &CosmologySimulationGridDimension[gridnum][1],
		    &CosmologySimulationGridDimension[gridnum][2]);
    if (sscanf(line, "CosmologySimulationGridLevel[%"ISYM"]", &gridnum) > 0)
      ret += sscanf(line, "CosmologySimulationGridLevel[%"ISYM"] = %"ISYM,
		    &gridnum, &CosmologySimulationGridLevel[gridnum]);
 
    ret += sscanf(line, "CosmologySimulationInitialFractionHII = %"FSYM,
		  &CosmologySimulationInitialFractionHII);
    ret += sscanf(line, "CosmologySimulationInitialFractionHeII = %"FSYM,
		  &CosmologySimulationInitialFractionHeII);
    ret += sscanf(line, "CosmologySimulationInitialFractionHeIII = %"FSYM,
		  &CosmologySimulationInitialFractionHeIII);
    ret += sscanf(line, "CosmologySimulationInitialFractionHM = %"FSYM,
		  &CosmologySimulationInitialFractionHM);
    ret += sscanf(line, "CosmologySimulationInitialFractionH2I = %"FSYM,
		  &CosmologySimulationInitialFractionH2I);
    ret += sscanf(line, "CosmologySimulationInitialFractionH2II = %"FSYM,
		  &CosmologySimulationInitialFractionH2II);
    ret += sscanf(line, "CosmologySimulationUseMetallicityField = %"ISYM,
		  &CosmologySimulationUseMetallicityField);
 
    ret += sscanf(line, "CosmologySimulationManuallySetParticleMassRatio = %"ISYM,
		  &CosmologySimulationManuallySetParticleMassRatio);
    ret += sscanf(line, "CosmologySimulationManualParticleMassRatio = %"FSYM,
		  &CosmologySimulationManualParticleMassRatio);

    ret += sscanf(line, "CosmologySimulationCalculatePositions = %"ISYM,
		  &CosmologySimulationCalculatePositions);

    // If the dummy char space was used, then make another
 
    if (dummy[0] != 0) {
      dummy = new char[MAX_LINE_LENGTH];
      ret++;
    }
 
    // if the line is suspicious, issue a warning
 
    if (ret == 0 && strstr(line, "=") && strstr(line, "CosmologySimulation") && 
	line[0] != '#' && MyProcessorNumber == ROOT_PROCESSOR)
      fprintf(stderr, "warning: the following parameter line was not interpreted:\n%s\n", line);
 
  }
 
  // More error checking
 
  if (CosmologySimulationDensityName == NULL &&
      (CosmologySimulationParticlePositionName == NULL &&
       !CosmologySimulationCalculatePositions)) {
    ENZO_FAIL("Missing initial data.\n");
  }
 
  if (CosmologySimulationDensityName != NULL && CellFlaggingMethod[0] != 2)
      fprintf(stderr, "CosmologySimulation: check CellFlaggingMethod.\n");
 
  if (CosmologySimulationDensityName == NULL && CellFlaggingMethod[0] != 4)
      fprintf(stderr, "CosmologySimulation: check CellFlaggingMethod.\n");
 
  if (CosmologySimulationNumberOfInitialGrids > MAX_INITIAL_GRIDS) {
    ENZO_FAIL("Too many InitialGrids! increase MAX_INITIAL_GRIDS\n");
  }
 
  if (CosmologySimulationDensityName == NULL && MultiSpecies+RadiativeCooling > 0) {
    fprintf(stderr, "warning: no density field; setting MultiSpecies/RadiativeCooling = 0\n");
    MultiSpecies = RadiativeCooling = 0;
  }

  if (CosmologySimulationParticleVelocityNames[0] != NULL &&
      !CosmologySimulationCalculatePositions) {
    ENZO_FAIL("CosmologySimulation: 1-component files only valid for use with "
	    "CosmologySimualtionCalculatePositions.\n");
  }
  // If temperature is left unset, set it assuming that T=550 K at z=200
 
  if (CosmologySimulationInitialTemperature == FLOAT_UNDEFINED)
    CosmologySimulationInitialTemperature = 550.0 *
      POW((1.0 + InitialRedshift)/(1.0 + 200), 2);
 
  // Generate the grids and set-up the hierarchy
 
  HierarchyEntry *GridsList[MAX_INITIAL_GRIDS];
  GridsList[0] = &TopGrid;
 
  if (MyProcessorNumber == ROOT_PROCESSOR)
    printf("Start loop creating initial grid hierarchy entries, from one\n");
 
  for (gridnum = 1; gridnum < CosmologySimulationNumberOfInitialGrids; gridnum++) {
 
    if (MyProcessorNumber == ROOT_PROCESSOR)
      printf("Create hierarchy entry for initial grid %"ISYM"\n", gridnum);
 
    // Create a spot in the hierarchy
 
    Subgrid    = new HierarchyEntry;
 
    // Find where to put this new grid
 
    int ParentGrid = INT_UNDEFINED;
 
    for (i = 0; i < gridnum; i++)
      if (CosmologySimulationGridLevel[i] ==
	  CosmologySimulationGridLevel[gridnum]-1)
	for (dim = 0; dim < MetaData.TopGridRank; dim++) {
	  if (CosmologySimulationGridLeftEdge[gridnum][dim] <
	      CosmologySimulationGridLeftEdge[i][dim]       ||
	      CosmologySimulationGridRightEdge[gridnum][dim] >
	      CosmologySimulationGridRightEdge[i][dim]       )
	    break;
	  ParentGrid = i;
	}
 
    if (ParentGrid == INT_UNDEFINED) {
      ENZO_VFAIL("Grid %"ISYM" has no valid parent.\n", gridnum)
    }
 
    // Insert this grid at the appropriate position in the subgrid chain
 
    GridsList[gridnum] = Subgrid;
    Subgrid->NextGridNextLevel = NULL;
    Subgrid->NextGridThisLevel = GridsList[ParentGrid]->NextGridNextLevel;
    Subgrid->ParentGrid        = GridsList[ParentGrid];
    GridsList[ParentGrid]->NextGridNextLevel = Subgrid;
 
    // Error check for consistency and add ghost zones to dimension
 
    for (dim = 0; dim < MetaData.TopGridRank; dim++) {
      FLOAT SubgridCellSize = (DomainRightEdge[dim] - DomainLeftEdge[dim])/
	FLOAT(MetaData.TopGridDims[dim]*
	      POW(FLOAT(RefineBy), CosmologySimulationGridLevel[gridnum]));

      if (debug1) {
	printf("  %"GSYM"\n", SubgridCellSize);
	printf("  %"GSYM"\n", POW(FLOAT(RefineBy), CosmologySimulationGridLevel[gridnum]));
	printf("  %"ISYM" %"ISYM"\n", MetaData.TopGridDims[dim], dim);
	printf("  %"GSYM" %"GSYM"\n", CosmologySimulationGridRightEdge[gridnum][dim],
	       CosmologySimulationGridLeftEdge[gridnum][dim]);
	printf("  %"ISYM"\n", nint((CosmologySimulationGridRightEdge[gridnum][dim] -
				    CosmologySimulationGridLeftEdge[gridnum][dim]   )
				   /SubgridCellSize));
      } // ENDIF debug1
 
      // Check if declared size matches left/right edges
 
      if (nint((CosmologySimulationGridRightEdge[gridnum][dim] -
		CosmologySimulationGridLeftEdge[gridnum][dim]   )
	       /SubgridCellSize) != CosmologySimulationGridDimension[gridnum][dim]) {
	fprintf(stderr, "Subgrid inconsistency: grid %"ISYM", dim %"ISYM"\n",
		gridnum, dim);
	fprintf(stderr, " subgrid: %"GOUTSYM" -> %"GOUTSYM", CellSize = %"GOUTSYM"\n",
	      CosmologySimulationGridLeftEdge[gridnum][dim],
	      CosmologySimulationGridRightEdge[gridnum][dim], SubgridCellSize);
	ENZO_FAIL("Subgrid Inconsistency!\n");
      }
 
      // Check if left/right edge fall on Parent cell boundary
 
      if (nint((CosmologySimulationGridLeftEdge[gridnum][dim] -
		CosmologySimulationGridLeftEdge[ParentGrid][dim])/
	       SubgridCellSize) % RefineBy != 0 ||
	  nint((CosmologySimulationGridRightEdge[gridnum][dim] -
		CosmologySimulationGridLeftEdge[ParentGrid][dim])/
	       SubgridCellSize) % RefineBy != 0 ) {
	fprintf(stderr, "Subgrid inconsistency: grid %"ISYM", dim %"ISYM"\n",
		gridnum, dim);
	fprintf(stderr, "left or right edges are not on parent cell edge.\n");
	ENZO_FAIL("Subgrid Inconsistency!\n");
      }
 
      // Add ghost zones
 
      CosmologySimulationGridDimension[gridnum][dim] += 2*DEFAULT_GHOST_ZONES;
 
    } // end of loop over dimensions
 
 
    // Create a new subgrid and initialize it
 
    Subgrid->GridData = new grid;
    Subgrid->GridData->InheritProperties(Subgrid->ParentGrid->GridData);
    Subgrid->GridData->PrepareGrid(MetaData.TopGridRank,
				   CosmologySimulationGridDimension[gridnum],
				   CosmologySimulationGridLeftEdge[gridnum],
				   CosmologySimulationGridRightEdge[gridnum],
				   0);
 
    // If subgrids are static, convert to static regions
 
    if (CosmologySimulationSubgridsAreStatic == TRUE) {
      for (region = 0; region < MAX_STATIC_REGIONS; region++)
	if (StaticRefineRegionLevel[region] == INT_UNDEFINED) {
	  StaticRefineRegionLevel[region] =
	    CosmologySimulationGridLevel[gridnum] - 1;
	  for (dim = 0; dim < MetaData.TopGridRank; dim++) {
	    StaticRefineRegionLeftEdge[region][dim] =
	      CosmologySimulationGridLeftEdge[gridnum][dim];
	    StaticRefineRegionRightEdge[region][dim] =
	      CosmologySimulationGridRightEdge[gridnum][dim];
	  }
	  for (dim = MetaData.TopGridRank; dim < MAX_DIMENSION; dim++) {
	    StaticRefineRegionLeftEdge[region][dim] = DomainLeftEdge[dim];
	    StaticRefineRegionRightEdge[region][dim] = DomainRightEdge[dim];
	  }
	  break;
	}
      if (region == MAX_STATIC_REGIONS) {
	ENZO_FAIL("Increase number of static refine regions\n");
      }
    }
 
    // Remove ghost zones from dim
 
    for (dim = 0; dim < MetaData.TopGridRank; dim++)
      CosmologySimulationGridDimension[gridnum][dim] -= 2*DEFAULT_GHOST_ZONES;
 
  } // end: loop over gridnums
 
 
  //---------------------------------------------------------------------------
 
 
  if (MyProcessorNumber == ROOT_PROCESSOR)
    printf("Start loop initializing generated grids, from zero\n");
 
  // Initialize the previously-generated grids
 
  for (gridnum = 0; gridnum < CosmologySimulationNumberOfInitialGrids; gridnum++) {
 
    if (MyProcessorNumber == ROOT_PROCESSOR)
      printf("RH: CosmologySimulation: Initializing grid %"ISYM"\n", gridnum);
 
    // If there is more than one grid, add the grid number to the name
 
    if (CosmologySimulationNumberOfInitialGrids > 1) {
 
      if (MyProcessorNumber == ROOT_PROCESSOR)
	printf("CosmologySimulation: Initializing grid %"ISYM"\n", gridnum);
 
      if (CosmologySimulationDensityName)
	sprintf(DensityName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationDensityName, gridnum);
      if (CosmologySimulationTotalEnergyName)
	sprintf(TotalEnergyName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationTotalEnergyName, gridnum);
      if (CosmologySimulationGasEnergyName)
	sprintf(GasEnergyName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationGasEnergyName, gridnum);
      for (dim = 0; dim < MetaData.TopGridRank; dim++)
	if (CosmologySimulationVelocityNames[dim])
	  sprintf(VelocityNames[dim] = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		  CosmologySimulationVelocityNames[dim], gridnum);
      if (CosmologySimulationParticlePositionName)
	sprintf(ParticlePositionName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationParticlePositionName, gridnum);
      if (CosmologySimulationParticleVelocityName)
	sprintf(ParticleVelocityName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationParticleVelocityName, gridnum);
      if (CosmologySimulationParticleMassName)
	sprintf(ParticleMassName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationParticleMassName, gridnum);
      if (CosmologySimulationParticleTypeName)
        sprintf(ParticleTypeName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
                CosmologySimulationParticleTypeName, gridnum);
      for (dim = 0; dim < MetaData.TopGridRank; dim++)
	if (CosmologySimulationParticleVelocityNames[dim])
	  sprintf(ParticleVelocityNames[dim] = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		  CosmologySimulationParticleVelocityNames[dim], gridnum);
	  
 
    } else {
 
      DensityName            = CosmologySimulationDensityName;
      TotalEnergyName        = CosmologySimulationTotalEnergyName;
      GasEnergyName          = CosmologySimulationGasEnergyName;
      for (dim = 0; dim < MetaData.TopGridRank; dim++)
	VelocityNames[dim]   = CosmologySimulationVelocityNames[dim];
      ParticlePositionName   = CosmologySimulationParticlePositionName;
      ParticleVelocityName   = CosmologySimulationParticleVelocityName;
      ParticleMassName       = CosmologySimulationParticleMassName;
      ParticleTypeName       = CosmologySimulationParticleTypeName;
      for (dim = 0; dim < MetaData.TopGridRank; dim++)
	ParticleVelocityNames[dim]   = CosmologySimulationParticleVelocityNames[dim];

    }
 
    // If there is a subgrid, use CosmologySimulationSubgridsAreStatic,
    // otherwise just set to false
 
    SubgridsAreStatic = (GridsList[gridnum]->NextGridNextLevel == NULL) ?
      FALSE : CosmologySimulationSubgridsAreStatic;
 
    // Initialize the grid by reading in (no) data
 
    int TotalRefinement = nint(POW(FLOAT(RefineBy),
				   CosmologySimulationGridLevel[gridnum]));
 
    if (MyProcessorNumber == ROOT_PROCESSOR)
      printf("Call CSIG for gridnum %"ISYM" with TR %"ISYM" and Dname %s\n", gridnum, TotalRefinement, DensityName);
 
    if (GridsList[gridnum]->GridData->NestedCosmologySimulationInitializeGrid(
			     gridnum,
			     CosmologySimulationOmegaBaryonNow,
			       CosmologySimulationOmegaCDMNow,
			       CosmologySimulationInitialTemperature,
			     DensityName, TotalEnergyName,
			       GasEnergyName, VelocityNames,
			       ParticlePositionName, ParticleVelocityName,
			       ParticleMassName, ParticleTypeName,
			     ParticleVelocityNames,
			     SubgridsAreStatic, TotalRefinement,
			     CosmologySimulationInitialFractionHII,
			     CosmologySimulationInitialFractionHeII,
			     CosmologySimulationInitialFractionHeIII,
			     CosmologySimulationInitialFractionHM,
			     CosmologySimulationInitialFractionH2I,
			     CosmologySimulationInitialFractionH2II,
			     CosmologySimulationUseMetallicityField,
			     MetaData.NumberOfParticles,
			     CosmologySimulationManuallySetParticleMassRatio,
			     CosmologySimulationManualParticleMassRatio,
			     CosmologySimulationCalculatePositions
						       ) == FAIL) {
      ENZO_FAIL("Error in grid->NestedCosmologySimulationInitializeGrid.\n");
    }
 
    // Set boundary conditions if necessary
 
  } // end loop over initial grids
 
  if (MyProcessorNumber == ROOT_PROCESSOR)
    printf("End of loop over initial grids\n");
 
  //---------------------------------------------------------------------------
 
 
  /* Convert minimum initial overdensity for refinement to mass
     (unless MinimumMass itself was actually set).
     Note: multiply MinimumMassForRefinement by the OmegaBaryonNow since the
     routine that uses this parameter only counts baryonic mass. */
 
  for (i = 0; i < MAX_FLAGGING_METHODS; i++) {
    if (MinimumMassForRefinement[i] == FLOAT_UNDEFINED) {
 
      MinimumMassForRefinement[i] = CosmologySimulationOmegaBaryonNow/
	                            OmegaMatterNow;
      if (CellFlaggingMethod[i] == 4)
	MinimumMassForRefinement[i] = CosmologySimulationOmegaCDMNow/
	                              OmegaMatterNow;
 
      MinimumMassForRefinement[i] *= MinimumOverDensityForRefinement[i];
      for (dim = 0; dim < MetaData.TopGridRank; dim++)
	MinimumMassForRefinement[i] *=
	  (DomainRightEdge[dim]-DomainLeftEdge[dim])/
	  float(MetaData.TopGridDims[dim]);
    }
  }
 
  // set up field names and units
 
  i = 0;
  DataLabel[i++] = DensName;
  DataLabel[i++] = TEName;
  if (DualEnergyFormalism)
    DataLabel[i++] = GEName;
  DataLabel[i++] = Vel1Name;
  DataLabel[i++] = Vel2Name;
  DataLabel[i++] = Vel3Name;
  if (MultiSpecies) {
    DataLabel[i++] = ElectronName;
    DataLabel[i++] = HIName;
    DataLabel[i++] = HIIName;
    DataLabel[i++] = HeIName;
    DataLabel[i++] = HeIIName;
    DataLabel[i++] = HeIIIName;
    if (MultiSpecies > 1) {
      DataLabel[i++] = HMName;
      DataLabel[i++] = H2IName;
      DataLabel[i++] = H2IIName;
    }
    if (MultiSpecies > 2) {
      DataLabel[i++] = DIName;
      DataLabel[i++] = DIIName;
      DataLabel[i++] = HDIName;
    }
  }
 
  if (CosmologySimulationUseMetallicityField) {
    DataLabel[i++] = MetalName;
    if(MultiMetals){
      DataLabel[i++] = ExtraNames[0];
      DataLabel[i++] = ExtraNames[1];
    }
  }
 
  if(STARMAKE_METHOD(COLORED_POP3_STAR)){
    DataLabel[i++] = ForbidName;
  }

  if (WritePotential)
    DataLabel[i++] = GPotName;
 
  for (j = 0; j < i; j++)
    DataUnits[j] = NULL;
 
  // Write parameters to parameter output file
 
  if (MyProcessorNumber == ROOT_PROCESSOR) {
    fprintf(Outfptr, "CosmologySimulationOmegaBaryonNow       = %"FSYM"\n",
	    CosmologySimulationOmegaBaryonNow);
    fprintf(Outfptr, "CosmologySimulationOmegaCDMNow          = %"FSYM"\n",
	    CosmologySimulationOmegaCDMNow);
    fprintf(Outfptr, "CosmologySimulationInitialTemperature   = %"FSYM"\n\n",
	    CosmologySimulationInitialTemperature);
 
    fprintf(Outfptr, "CosmologySimulationDensityName          = %s\n",
	    CosmologySimulationDensityName);
    if (CosmologySimulationTotalEnergyName)
    fprintf(Outfptr, "CosmologySimulationTotalEnergyName      = %s\n",
	    CosmologySimulationTotalEnergyName);
    if (CosmologySimulationGasEnergyName)
    fprintf(Outfptr, "CosmologySimulationGasEnergyName        = %s\n",
	    CosmologySimulationGasEnergyName);
    fprintf(Outfptr, "CosmologySimulationVelocity1Name        = %s\n",
	    CosmologySimulationVelocityNames[0]);
    fprintf(Outfptr, "CosmologySimulationVelocity2Name        = %s\n",
	    CosmologySimulationVelocityNames[1]);
    fprintf(Outfptr, "CosmologySimulationVelocity3Name        = %s\n",
	    CosmologySimulationVelocityNames[2]);
    if (CosmologySimulationParticlePositionName)
    fprintf(Outfptr, "CosmologySimulationParticlePositionName = %s\n",
	    CosmologySimulationParticlePositionName);
    if (CosmologySimulationParticleVelocityName)
    fprintf(Outfptr, "CosmologySimulationParticleVelocityName = %s\n",
	    CosmologySimulationParticleVelocityName);
    if (CosmologySimulationParticleMassName)
    fprintf(Outfptr, "CosmologySimulationParticleMassName     = %s\n",
	    CosmologySimulationParticleMassName);
    if (CosmologySimulationParticleTypeName)
    fprintf(Outfptr, "CosmologySimulationParticleTypeName     = %s\n\n",
            CosmologySimulationParticleTypeName);
    fprintf(Outfptr, "CosmologySimulationParticleVelocity1Name = %s\n",
	    CosmologySimulationVelocityNames[0]);
    fprintf(Outfptr, "CosmologySimulationParticleVelocity2Name = %s\n",
	    CosmologySimulationVelocityNames[1]);
    fprintf(Outfptr, "CosmologySimulationParticleVelocity3Name = %s\n",
	    CosmologySimulationVelocityNames[2]);
 
    fprintf(Outfptr, "CosmologySimulationNumberOfInitialGrids = %"ISYM"\n",
	    CosmologySimulationNumberOfInitialGrids);
    fprintf(Outfptr, "CosmologySimulationSubgridsAreStatic    = %"ISYM"\n",
	    CosmologySimulationSubgridsAreStatic);
    fprintf(Outfptr, "CosmologySimulationCalculatePositions   = %"ISYM"\n",
	    CosmologySimulationCalculatePositions);
    
 
    for (gridnum = 1; gridnum < CosmologySimulationNumberOfInitialGrids;
	 gridnum++) {
      fprintf(Outfptr, "CosmologySimulationGridLeftEdge[%"ISYM"]     = ", gridnum);
      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
			CosmologySimulationGridLeftEdge[gridnum]);
      fprintf(Outfptr, "CosmologySimulationGridRightEdge[%"ISYM"]    = ", gridnum);
      WriteListOfFloats(Outfptr, MetaData.TopGridRank,
			CosmologySimulationGridRightEdge[gridnum]);
      fprintf(Outfptr, "CosmologySimulationGridDimension[%"ISYM"]    = ", gridnum);
      WriteListOfInts(Outfptr, MetaData.TopGridRank,
		      CosmologySimulationGridDimension[gridnum]);
    }
 
    fprintf(Outfptr, "\n");
 
    fprintf(Outfptr, "CosmologySimulationInitialFractionHII   = %"GSYM"\n",
	    CosmologySimulationInitialFractionHII);
    fprintf(Outfptr, "CosmologySimulationInitialFractionHeII  = %"GSYM"\n",
	    CosmologySimulationInitialFractionHeII);
    fprintf(Outfptr, "CosmologySimulationInitialFractionHeIII = %"GSYM"\n",
	    CosmologySimulationInitialFractionHeIII);
    fprintf(Outfptr, "CosmologySimulationInitialFractionHM    = %"GSYM"\n",
	    CosmologySimulationInitialFractionHM);
    fprintf(Outfptr, "CosmologySimulationInitialFractionH2I   = %"GSYM"\n",
	    CosmologySimulationInitialFractionH2I);
    fprintf(Outfptr, "CosmologySimulationInitialFractionH2II  = %"GSYM"\n",
	    CosmologySimulationInitialFractionH2II);
    fprintf(Outfptr, "CosmologySimulationUseMetallicityField  = %"ISYM"\n\n",
	    CosmologySimulationUseMetallicityField);
  }
 
  // Clean up
 
  delete [] dummy;
 
  return SUCCESS;
}
 
 
 
 
void NestedRecursivelySetParticleCount(HierarchyEntry *GridPoint, PINT *Count);
 
 
// Re-call the initializer on level zero grids.
// Used in case of ParallelRootGridIO.
 
 
int NestedCosmologySimulationReInitialize(HierarchyEntry *TopGrid,
				          TopGridData &MetaData)
{
 
 
  int dim, gridnum = 0;
 
  HierarchyEntry *CurrentGrid;
  HierarchyEntry *Temp;
 
 
  char *DensityName = NULL, *TotalEnergyName = NULL, *GasEnergyName = NULL,
    *ParticlePositionName = NULL, *ParticleVelocityName = NULL,
    *ParticleMassName = NULL, *VelocityNames[MAX_DIMENSION],
    *ParticleTypeName = NULL, *ParticleVelocityNames[MAX_DIMENSION];;
 
 for (dim = 0; dim < MAX_DIMENSION; dim++) {
    ParticleVelocityNames[dim] = NULL;
    VelocityNames[dim] = NULL;
  }
 
  CurrentGrid = TopGrid;
 
  for (gridnum = 0; gridnum < CosmologySimulationNumberOfInitialGrids; gridnum++) {
 
    if (MyProcessorNumber == ROOT_PROCESSOR)
      printf("NestedCosmologySimulation: ReInitializing grid %"ISYM"\n", gridnum);
 
  // If there is more than one grid, add the grid number to the name
 
  if (CosmologySimulationNumberOfInitialGrids > 1) {
 
    if (CosmologySimulationDensityName)
      sprintf(DensityName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationDensityName, gridnum);
    if (CosmologySimulationTotalEnergyName)
      sprintf(TotalEnergyName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationTotalEnergyName, gridnum);
    if (CosmologySimulationGasEnergyName)
      sprintf(GasEnergyName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationGasEnergyName, gridnum);
 
    for (dim = 0; dim < MAX_DIMENSION; dim++)
      if (CosmologySimulationVelocityNames[dim])
	sprintf(VelocityNames[dim] = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationVelocityNames[dim], gridnum);
 
    if (CosmologySimulationParticlePositionName)
      sprintf(ParticlePositionName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationParticlePositionName, gridnum);
    if (CosmologySimulationParticleVelocityName)
      sprintf(ParticleVelocityName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationParticleVelocityName, gridnum);
    if (CosmologySimulationParticleMassName)
      sprintf(ParticleMassName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
	      CosmologySimulationParticleMassName, gridnum);
    if (CosmologySimulationParticleTypeName)
      sprintf(ParticleTypeName = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
              CosmologySimulationParticleTypeName, gridnum);
    for (dim = 0; dim < MAX_DIMENSION; dim++)
      if (CosmologySimulationParticleVelocityNames[dim])
	sprintf(ParticleVelocityNames[dim] = new char[MAX_LINE_LENGTH], "%s.%1"ISYM,
		CosmologySimulationParticleVelocityNames[dim], gridnum);
 
  } else {
 
    DensityName            = CosmologySimulationDensityName;
    TotalEnergyName        = CosmologySimulationTotalEnergyName;
    GasEnergyName          = CosmologySimulationGasEnergyName;
    for (dim = 0; dim < MAX_DIMENSION; dim++)
      VelocityNames[dim]   = CosmologySimulationVelocityNames[dim];
    ParticlePositionName   = CosmologySimulationParticlePositionName;
    ParticleVelocityName   = CosmologySimulationParticleVelocityName;
    ParticleMassName       = CosmologySimulationParticleMassName;
    ParticleTypeName       = CosmologySimulationParticleTypeName;
    for (dim = 0; dim < MAX_DIMENSION; dim++)
      ParticleVelocityNames[dim]   = CosmologySimulationParticleVelocityNames[dim];
 
  }
 
  // If there is a subgrid, use CosmologySimulationSubgridsAreStatic,
  // otherwise just set to false
 
  int SubgridsAreStatic = (CurrentGrid->NextGridNextLevel == NULL) ?
    FALSE : CosmologySimulationSubgridsAreStatic;
 
  // Call grid initializer.  Use TotalRefinement = -1 to flag real read
 
  int TotalRefinement = -1;
 
  // Loop over level zero grid
 
  Temp = CurrentGrid;
 
  while (Temp != NULL) {
 
    if (Temp->GridData->NestedCosmologySimulationInitializeGrid(gridnum,
			     CosmologySimulationOmegaBaryonNow,
			       CosmologySimulationOmegaCDMNow,
			       CosmologySimulationInitialTemperature,
			     DensityName, TotalEnergyName,
			       GasEnergyName, VelocityNames,
			       ParticlePositionName, ParticleVelocityName,
			       ParticleMassName, ParticleTypeName,
			     ParticleVelocityNames,
			     SubgridsAreStatic, TotalRefinement,
			     CosmologySimulationInitialFractionHII,
			     CosmologySimulationInitialFractionHeII,
			     CosmologySimulationInitialFractionHeIII,
			     CosmologySimulationInitialFractionHM,
			     CosmologySimulationInitialFractionH2I,
			     CosmologySimulationInitialFractionH2II,
			     CosmologySimulationUseMetallicityField,
			     MetaData.NumberOfParticles,
			     CosmologySimulationManuallySetParticleMassRatio,
			     CosmologySimulationManualParticleMassRatio,
			     CosmologySimulationCalculatePositions
						       ) == FAIL) {
      ENZO_FAIL("Error in grid->NestedCosmologySimulationInitializeGrid.\n");
    }
 
    Temp = Temp->NextGridThisLevel;
  }
 
  CurrentGrid = CurrentGrid->NextGridNextLevel;
 
  } // end loop over initial grid levels
 
  // Create tracer particles
 
  PINT DummyNumberOfParticles = 0;
 
  Temp = TopGrid;
 
  while (Temp != NULL) {
    if (Temp->GridData->TracerParticleCreateParticles(
                                     TracerParticleCreationLeftEdge,
                                     TracerParticleCreationRightEdge,
                                     TracerParticleCreationSpacing,
                                     DummyNumberOfParticles) == FAIL) {
      ENZO_FAIL("Error in grid->TracerParticleCreateParticles\n");
    }
 
    Temp = Temp->NextGridThisLevel;
  }
 
 
  // Get the global particle count
 
  int LocalNumberOfParticles;
 
  CurrentGrid = TopGrid;
 
  while (CurrentGrid != NULL) {
 
    Temp = CurrentGrid;
 
    while (Temp != NULL) {
 
      LocalNumberOfParticles = 0;
      LocalNumberOfParticles = Temp->GridData->ReturnNumberOfParticles();
      // printf("OldLocalParticleCount: %"ISYM"\n", LocalNumberOfParticles );
 
#ifdef USE_MPI
      CommunicationAllReduceValues(&LocalNumberOfParticles, 1, MPI_SUM);
#endif /* USE_MPI */
      Temp->GridData->SetNumberOfParticles(LocalNumberOfParticles);
 
      LocalNumberOfParticles = Temp->GridData->ReturnNumberOfParticles();
      // printf("NewLocalParticleCount: %"ISYM"\n", LocalNumberOfParticles );
 
      Temp = Temp->NextGridThisLevel;
    }
 
    CurrentGrid = CurrentGrid->NextGridNextLevel;
 
  }
 
  // Loop over grids and set particle ID number
 
  Temp = TopGrid;
  PINT ParticleCount = 0;
 
  NestedRecursivelySetParticleCount(Temp, &ParticleCount);
 
  if (debug)
    printf("FinalParticleCount = %"PISYM"\n", ParticleCount);
 
  MetaData.NumberOfParticles = ParticleCount;

  /* Now set the grids' parents correctly instead of having all of the
     siblings have the same parent.  Before this left the work for
     RebuildHierarchy, which will send all of the regions and
     particles back to one processor! */

  bool match;
  int level, Rank, TempDims[MAX_DIMENSION];
  FLOAT GridCenter[MAX_DIMENSION];
  LevelHierarchyEntry *LevelArray[MAX_DEPTH_OF_HIERARCHY];
  LevelHierarchyEntry *CArray, *PArray;
  HierarchyEntry *Current, *Parent;  // For convenience
  FLOAT LeftEdge[MAX_DIMENSION], RightEdge[MAX_DIMENSION];
  FLOAT LeftParent[MAX_DIMENSION], RightParent[MAX_DIMENSION];

  // Level arrays are convenient
  for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++)
    LevelArray[level] = NULL;
  AddLevel(LevelArray, TopGrid, 0);

  /* Now that we have the pointers stored in LevelArray, let's reset
     all of the hierarchy pointers (except for the level-0
     NextGridThisLevel) */

  for (level = 0; level < MAX_DEPTH_OF_HIERARCHY; level++) {
    for (CArray = LevelArray[level]; CArray; CArray = CArray->NextGridThisLevel) {
      if (level > 0)
	CArray->GridHierarchyEntry->NextGridThisLevel = NULL;
      CArray->GridHierarchyEntry->NextGridNextLevel = NULL;
      CArray->GridHierarchyEntry->ParentGrid = NULL;
    } // ENDFOR grids
  } // ENDFOR level

  // Loop over all child grids, looking for parents
  for (level = 1; level < MAX_DEPTH_OF_HIERARCHY; level++) {
    for (CArray = LevelArray[level]; CArray; CArray = CArray->NextGridThisLevel) {

      Current = CArray->GridHierarchyEntry;
      Current->GridData->ReturnGridInfo(&Rank, TempDims, LeftEdge, RightEdge);
      for (dim = 0; dim < Rank; dim++)
	GridCenter[dim] = 0.5*(LeftEdge[dim] + RightEdge[dim]);

      // Look for parents in level-1
      for (PArray = LevelArray[level-1]; PArray; 
	   PArray = PArray->NextGridThisLevel) {

	Parent = PArray->GridHierarchyEntry;
	Parent->GridData->ReturnGridInfo(&Rank, TempDims, LeftParent, RightParent);

	// Grid won't necessarily be fully contained in initial
	// "parent".  Here we just look for overlap if the partitioned
	// grids weren't created to be fully contained within a
	// partition of the level-1 grid.

	match = true;
	for (dim = 0; dim < Rank; dim++)
	  match &= (GridCenter[dim] > LeftParent[dim]) &&
	    (GridCenter[dim] < RightParent[dim]);

	if (match) {
	  Current->ParentGrid = Parent;
	  Current->NextGridThisLevel = Parent->NextGridNextLevel;
	  Parent->NextGridNextLevel = Current;
	  break;
	} // ENDIF match

      } // ENDFOR parent level grids
    } // ENDFOR current level grids
  } // ENDFOR level

  return SUCCESS;
}
 
 
 
 
void NestedRecursivelySetParticleCount(HierarchyEntry *GridPoint, PINT *Count)
{
  // Add Count to the particle id's on this grid (which start from zero
  // since we are doing a parallel root grid i/o)
 
  GridPoint->GridData->AddToParticleNumber(Count);
 
  // Recursively apply this to siblings and children
 
  if (GridPoint->NextGridThisLevel != NULL)
    NestedRecursivelySetParticleCount(GridPoint->NextGridThisLevel, Count);
 
  if (GridPoint->NextGridNextLevel != NULL)

    NestedRecursivelySetParticleCount(GridPoint->NextGridNextLevel, Count);
 
  CommunicationBroadcastValue(Count, ROOT_PROCESSOR);
  return;
}
